###### Figure A12 Script #####
rm(list=ls());gc();gc();gc();gc();gc();gc();gc();gc()


library(PanelMatch) #using PanelMatch 3.0.0
library(dplyr)


here::i_am("Scripts/FigA12Script.R")

library(here)

kl_e2<-readRDS(here("Data", "FigureA12_data.rds"))
kl_e2$time<-NA;
kl_e2$time[which(kl_e2$yearsess2==1999)]<-1
kl_e2$time[which(kl_e2$yearsess2==2001)]<-2
kl_e2$time[which(kl_e2$yearsess2==2003)]<-3
kl_e2$time[which(kl_e2$yearsess2==2005)]<-4
kl_e2$time[which(kl_e2$yearsess2==2007)]<-5
kl_e2$time[which(kl_e2$yearsess2==2009)]<-6
kl_e2$time[which(kl_e2$yearsess2==2011)]<-7
kl_e2$time[which(kl_e2$yearsess2==2013)]<-8
kl_e2$time[which(kl_e2$yearsess2==2015)]<-9
kl_e2$time[which(kl_e2$yearsess2==2017)]<-10
kl_e2$time<-as.integer(kl_e2$time)
kl_e2<-kl_e2[order(kl_e2$CandIdYears),]

kl_e2$party2<-as.numeric(grepl("democrat",kl_e2$party)==T) #n.b. sole independent caucusing with reps (at least second time, first to be confirmed)

kl_e2$ky<-as.numeric(kl_e2$sab=="KY")
kl_e2$mi<-as.numeric(kl_e2$sab=="MI")
kl_e2$oh<-as.numeric(kl_e2$sab=="OH")
kl_e2$wv<-as.numeric(kl_e2$sab=="WV")

kl_e2$legislative_seniority<-as.numeric(kl_e2$incumbent_legislature)

kl_e2$sen<-as.numeric(kl_e2$sen)
kl_e2$leader<-as.numeric(kl_e2$leader)
kl_e2$chairs_<-as.numeric(kl_e2$chairs_)
kl_e2$CandIdSector<-as.integer(kl_e2$CandIdSector)

kl_e2$exc0<-as.numeric(kl_e2$sen==1 & kl_e2$yearsess2 ==2019 & kl_e2$sab=="WV" | kl_e2$sen==1 & kl_e2$yearsess2 ==2019 & kl_e2$sab=="OH" | kl_e2$sen==1 & kl_e2$yearsess2 ==2019 & kl_e2$sab=="KY" | kl_e2$cleaned_candid==262471 |kl_e2$filerid2=="NA")

kl_e2$sector_1<-as.numeric(kl_e2$sector=="agriculture")
kl_e2$sector_2<-as.numeric(kl_e2$sector=="business")
kl_e2$sector_3<-as.numeric(kl_e2$sector=="construction")
kl_e2$sector_4<-as.numeric(kl_e2$sector=="education")
kl_e2$sector_5<-as.numeric(kl_e2$sector=="energy")
kl_e2$sector_6<-as.numeric(kl_e2$sector=="finance")
kl_e2$sector_7<-as.numeric(kl_e2$sector=="health")
kl_e2$sector_8<-as.numeric(kl_e2$sector=="law")
kl_e2$sector_9<-as.numeric(kl_e2$sector=="transportation")

kl_e3<-as.data.frame(kl_e2,stringsAsFactors=F)

panel.data<-PanelData(panel.data=subset(kl_e3,limit1!=1&exc0==0), time.id = "time", unit.id = "CandIdSector", treatment = "chairs_",outcome="logfriamount_")

DiD.matches_chairs <- PanelMatch(lag = 1, 
                                  refinement.method = "CBPS.weight", 
                                 panel.data = panel.data , match.missing = T, 
                                 covs.formula =
                                   ~ I(lag(logfriamount_,1))+
                                   I(lag(mi,1))+ I(lag(oh,1))+I(lag(wv,1))+ I(lag(ky,1))+
                                   I(lag(sen,1))+  I(lag(nonelses,1))+I(lag(party2,1))+I(lag(incmaj,1)) +  I(lag(legislative_seniority,1))+  I(lag(leader,1))
                                 +   I(lag(sector_2,1))+ I(lag(sector_3,1))+ I(lag(sector_4,1))+ I(lag(sector_5,1))+
                                   I(lag(sector_6,1))+ I(lag(sector_7,1))+ I(lag(sector_8,1))+ I(lag(sector_9,1)),
                                 size.match = 20,
                                 qoi = "att",
                                 use.diagonal.variance.matrix = T,
                                  lead = 0, forbid.treatment.reversal = F)

#Estimation can take some time (generally less than an hour, but more than 10 minutes)
didresults2 <- PanelEstimate(confidence.level=c(0.9),sets = DiD.matches_chairs, panel.data = panel.data,number.iterations = 20000,df.adjustment = T)
summary(didresults2)  
pdf(file=here("Results", "FigA12.pdf"),width = 10,height = 6)
par(mfrow=c(1,2))
plot(didresults2,ylab="Estimated Effect of Treatment on Treated",main="(1) Chair Treatment for Sector j")


cb3<-get_covariate_balance(DiD.matches_chairs,panel.data=panel.data,covariates=c("incmaj","party2",
                                                                           "nonelses","legislative_seniority","leader","sen"))

plot(cb3[[1]]$att[,1],type="l",ylim=c(-2,2),ylab="SD",xaxt="n",xlab="Time",main="(2) Chair Treatment for Sector j -\n Covariate Balance",)
try(lines(cb3[[1]]$att[,2],col="red"))
try(lines(cb3[[1]]$att[,3],col="blue"))
try(lines(cb3[[1]]$att[,4],col="yellow"))
try(lines(cb3[[1]]$att[,5],col="green"))
try(lines(cb3[[1]]$att[,6],col="orange"))
axis(side=1,at=c(1,2),c("t-1","t-0"))
abline(h=0,lty=2)
legend("topleft",col=c("black","red","blue","yellow","green","orange"),lty=1,cex=0.9,legend=c("Majority Party","Party","Non-Election Cohort","Legislative Seniority","Leader","Senate"))
dev.off()

#This estimation involves some simulation, so results might be slightly different than what is 
#in the paper. To reproduce exactly the paper's version of these images, run the following code:
#(please note that commands in current version of PanelMatch, 3.0.0., require different commands to format and present data
#than the version of PanelMatch used for the original analyses, PanelMatch 2.0.0)

rm(list=ls());gc();gc();gc();gc();gc()

load(here("Data", "FigureA12.RData"))

#code producing objects DiD.matches_chairs and didresults2 based on object kl_e3
# DiD.matches_chairs <- PanelMatch(lag = 1, time.id = "time", unit.id = "CandIdSector", 
#                                  treatment = "chairs_", refinement.method = "CBPS.weight", 
#                                  data = subset(kl_e3,limit1!=1&exc0==0), match.missing = T, 
#                                  covs.formula =
#                                    ~ I(lag(logfriamount_,1))+
#                                    I(lag(mi,1))+ I(lag(oh,1))+I(lag(wv,1))+ I(lag(ky,1))+
#                                    I(lag(sen,1))+  I(lag(nonelses,1))+I(lag(party2,1))+I(lag(incmaj,1)) +  I(lag(legislative_seniority,1))+  I(lag(leader,1))
#                                  
#                                  +   I(lag(sector_2,1))+ I(lag(sector_3,1))+ I(lag(sector_4,1))+ I(lag(sector_5,1))+
#                                    I(lag(sector_6,1))+ I(lag(sector_7,1))+ I(lag(sector_8,1))+ I(lag(sector_9,1)),
#                                  size.match = 20,
#                                  qoi = "att",
#                                  use.diagonal.variance.matrix = T,
#                                  outcome.var = "logfriamount_", lead = 0, forbid.treatment.reversal = F)

#didresults2 <- PanelEstimate(confidence.level=c(0.9),sets = DiD.matches_chairs, data = subset(kl_e3,limit1!=1&exc0==0),number.iterations = 20000,df.adjustment = T)

summary(didresults2)  
#(exp(0.57)-1) * 100 = 76.82671
pdf(here("Results","FigA12Exact.pdf"),width = 10,height = 6)
par(mfrow=c(1,2))
plot(didresults2,ylab="Estimated Effect of Treatment on Treated",main="(1) Chair Treatment for Sector j")

panel.data<-PanelData(panel.data=subset(kl_e3,limit1!=1&exc0==0), time.id = "time", unit.id = "CandIdSector", treatment = "chairs_",outcome="logfriamount_")

cb3<-get_covariate_balance(DiD.matches_chairs,panel.data=panel.data,covariates=c("incmaj","party2",
                                                                                 "nonelses","legislative_seniority","leader","sen"))

plot(cb3[[1]]$att[,1],type="l",ylim=c(-2,2),ylab="SD",xaxt="n",xlab="Time",main="(2) Chair Treatment for Sector j -\n Covariate Balance",)
try(lines(cb3[[1]]$att[,2],col="red"))
try(lines(cb3[[1]]$att[,3],col="blue"))
try(lines(cb3[[1]]$att[,4],col="yellow"))
try(lines(cb3[[1]]$att[,5],col="green"))
try(lines(cb3[[1]]$att[,6],col="orange"))
axis(side=1,at=c(1,2),c("t-1","t-0"))
abline(h=0,lty=2)
legend("topleft",col=c("black","red","blue","yellow","green","orange"),lty=1,cex=0.9,legend=c("Majority Party","Party","Non-Election Cohort","Legislative Seniority","Leader","Senate"))

dev.off()

